Atelier R Cartographie

L’écosystème spatial sur R

Introduction à sf

  • Site web de sf: Simple Features for R

  • sf pour Simple Features

  • Sortie en octobre 2016

  • A pour but de rassembler les fonctionnalités d’anciens packages (sp, rgeos and rgdal) en un seul

  • Facilite la manipulation de données spatiales, avec des objets simples.

  • Tidy data: compatible avec la syntaxe pipe %>% et les opérateurs du tidyverse.

  • Principal auteur et mainteneur : Edzer Pebesma (également auteur du package sp)


la structure des objets sf :

format sf

Importer / exporter des données

library(sf)
## Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
mtq <- read_sf("data/mtq/martinique.shp")
mtq <- st_read("data/mtq/martinique.shp")
## Reading layer `martinique' from data source `/home/comeetie/Projets/quantilille/lecture/data/mtq/martinique.shp' using driver `ESRI Shapefile'
## Simple feature collection with 34 features and 23 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 690574.4 ymin: 1592426 xmax: 736126.5 ymax: 1645660
## Projected CRS: WGS 84 / UTM zone 20N
write_sf(mtq,"data/mtq/martinique.gpkg",delete_layer = TRUE)
st_write(mtq,"data/mtq/martinique.gpkg",delete_layer = TRUE)
## Deleting layer `martinique' using driver `GPKG'
## Writing layer `martinique' to data source `data/mtq/martinique.gpkg' using driver `GPKG'
## Writing 34 features with 23 fields and geometry type Polygon.

Projection

Obtenir la projection en utilisant st_crs() (code epsg) et la modifier en utilisant st_transform().

st_crs(mtq)
## Coordinate Reference System:
##   User input: WGS 84 / UTM zone 20N 
##   wkt:
## PROJCRS["WGS 84 / UTM zone 20N",
##     BASEGEOGCRS["WGS 84",
##         DATUM["World Geodetic System 1984",
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4326]],
##     CONVERSION["UTM zone 20N",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["Degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",-63,
##             ANGLEUNIT["Degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     ID["EPSG",32620]]
mtq_4326 <- mtq %>% st_transform(4326)

Afficher les données

Affichage par défaut :

plot(mtq)
## Warning: plotting the first 10 out of 23 attributes; use max.plot = 23 to plot
## all

En ne gardant que la géométrie :

plot(st_geometry(mtq))

Extraire les centroïdes

mtq_c <- st_centroid(mtq)
## Warning in st_centroid.sf(mtq): st_centroid assumes attributes are constant over
## geometries of x
plot(st_geometry(mtq))
plot(st_geometry(mtq_c), add=TRUE, cex=1.2, col="red", pch=20)

Matrice de distance

mat <- st_distance(x=mtq_c,y=mtq_c)
mat[1:5,1:5]
## Units: [m]
##           [,1]     [,2]      [,3]      [,4]      [,5]
## [1,]     0.000 35297.56  3091.501 12131.617 17136.310
## [2,] 35297.557     0.00 38332.602 25518.913 18605.249
## [3,]  3091.501 38332.60     0.000 15094.702 20226.198
## [4,] 12131.617 25518.91 15094.702     0.000  7177.011
## [5,] 17136.310 18605.25 20226.198  7177.011     0.000

Agrégation de polygones

Union simple :

mtq_u <- st_union(mtq)
plot(st_geometry(mtq), col="lightblue")
plot(st_geometry(mtq_u), add=T, lwd=2, border = "red")

A partir d’une variable de regroupement :

library(dplyr)
## 
## Attachement du package : 'dplyr'
## Les objets suivants sont masqués depuis 'package:stats':
## 
##     filter, lag
## Les objets suivants sont masqués depuis 'package:base':
## 
##     intersect, setdiff, setequal, union
mtq_u2 <- mtq %>% 
  group_by(STATUT) %>% 
  summarize(P13_POP = sum(P13_POP))
plot(st_geometry(mtq), col="lightblue")
plot(st_geometry(mtq_u2), add=T, lwd=2, border = "red", col=NA)

Zone tampon

mtq_b <- st_buffer(x = mtq_u, dist = 5000)
plot(st_geometry(mtq), col="lightblue")
plot(st_geometry(mtq_u), add=T, lwd=2)
plot(st_geometry(mtq_b), add=T, lwd=2, border = "red")

Intersection de polygones

# create a polygon
m <- rbind(c(700015,1624212), c(700015,1641586), c(719127,1641586), 
           c(719127,1624212), c(700015,1624212))
p <- st_sf(st_sfc(st_polygon(list(m))), crs = st_crs(mtq))
plot(st_geometry(mtq))
plot(p, border="red", lwd=2, add=T)

st_intersection() extrait la partie de mtq qui s’intersecte avec le polygone créé.

mtq_z <- st_intersection(x = mtq, y = p)
plot(st_geometry(mtq))
plot(st_geometry(mtq_z), col="red", border="green", add=T)

Compter les points dans des polygones

st_sample() crée des points aléatoires sur la carte de Paris.

pts <- st_sample(x = mtq, size = 50)
plot(st_geometry(mtq))
plot(pts, pch = 20, col = "red", add=TRUE, cex = 1)

st_interects() crée une liste de points dans chaque polygone.

inter <- st_intersects(mtq, pts)
mtq$nbpts <- sapply(X = inter, FUN = length)
plot(st_geometry(mtq))
# display munucipalities that intersect at least 2 point
plot(st_geometry(mtq[mtq$nbpts>2,]), col = "grey", add=TRUE)
plot(pts, pch = 20, col = "red", add=TRUE, cex = 1)

Autres packages

CRAN task views permet d’avoir des informations sur les packages du CRAN pertinents pour des tâches reliées à certains sujets.

CRAN Task View: Analysis of Spatial Data:

  • Classes for spatial data
  • Handling spatial data
  • Reading and writing spatial data
  • Visualisation
  • Point pattern analysis
  • Geostatistics
  • Disease mapping and areal data analysis
  • Spatial regression
  • Ecological analysis

Préparer / récupérer des données

library(sf)
library(dplyr)
# Import de la couche géographique (iris de Paris)
#iris.75 <- st_read(dsn = "../data/iris_75.shp",              stringsAsFactors = FALSE, quiet=TRUE)
iris.75 <- readRDS("../data/iris_75.RDS") %>% 
  st_set_crs(2154)
head(iris.75)
## Simple feature collection with 6 features and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 648979.6 ymin: 6861839 xmax: 654353.1 ymax: 6866159
## Projected CRS: RGF93 / Lambert-93
##     CODE_IRIS INSEE_COM                       geometry
## 78  751197316     75119 MULTIPOLYGON (((653970.6 68...
## 83  751176716     75117 MULTIPOLYGON (((649189.3 68...
## 177 751103703     75110 MULTIPOLYGON (((652767.6 68...
## 436 751187104     75118 MULTIPOLYGON (((652827.6 68...
## 468 751114314     75111 MULTIPOLYGON (((654272.9 68...
## 550 751103707     75110 MULTIPOLYGON (((652960.8 68...

Conversion csv vers sf

# Import du dataset  
accidents.2019.paris <- readRDS("../data/accidents2019_paris.RDS")
# Transformation en objet sf
accidents.2019.paris <- st_as_sf(accidents.2019.paris,
                                coords = c("long", "lat"),
                                crs = 4326, agr = "constant") %>% 
  st_transform(2154)
plot(st_geometry(accidents.2019.paris))

inter <- st_intersects(iris.75, accidents.2019.paris)
inter_blessgravtues <- st_intersects(iris.75, accidents.2019.paris
                              %>% filter(grav%in%c(2,3)))
iris.75$nbacc <- sapply(X = inter, FUN = length)
iris.75$nbacc_blessgravtues <- sapply(X = inter_blessgravtues, FUN = length)

#Il manque 24 accidents
nrow(accidents.2019.paris)-sum(iris.75$nbacc)
## [1] 24
plot(st_geometry(iris.75))
plot(accidents.2019.paris %>% filter(grav%in%c(1,4)) ,
     pch = 20, col = "darkgreen", add=TRUE, cex = 0.5)
## Warning in plot.sf(accidents.2019.paris %>% filter(grav %in% c(1, 4)), pch =
## 20, : ignoring all but the first attribute
plot(accidents.2019.paris %>% filter(grav%in%c(3)) ,
     pch = 20, col = "orange", add=TRUE, cex = 0.5)
## Warning in plot.sf(accidents.2019.paris %>% filter(grav %in% c(3)), pch = 20, :
## ignoring all but the first attribute
plot(accidents.2019.paris %>% filter(grav==2),
     pch = 20, col = "red", add=TRUE, cex = 1)
## Warning in plot.sf(accidents.2019.paris %>% filter(grav == 2), pch = 20, :
## ignoring all but the first attribute

Utiliser osmdata

osmdata permet d’extraire des éléments de la base de données gratuite et open-source OpenStreetMap.

library(osmdata)
# Récupérer les routes principales grâce à osm
bb      <- iris.75 %>% st_transform(4326) %>% st_bbox()
q       <- opq(bbox = bb,timeout = 180)
qm      <- add_osm_feature (q, key = 'highway',value = 'motorway', value_exact = FALSE)
qt      <- add_osm_feature (q, key = 'highway',value = 'trunk', value_exact = FALSE)
qp      <- add_osm_feature (q, key = 'highway',value = 'primary', value_exact = FALSE)

motorway<- osmdata_sf(qm)
trunk   <- osmdata_sf(qt)
primary <- osmdata_sf(qp)

roads    <- c(primary,trunk,motorway)$osm_lines %>% st_transform(st_crs(iris.75))
roads.geom = st_intersection(st_geometry(roads),iris.75) 

# Récupérer le shape de la seine 
qr <- q %>% 
  add_osm_feature (key = 'waterway') %>% 
  add_osm_feature(key = "name:fr", value = "La Seine")
river <- osmdata_sf(qr)

river.geom <- c(st_geometry(river$osm_lines),st_geometry(river$osm_multilines)) %>% st_transform(st_crs(iris.75))
#river.geom <- st_geometry(river$osm_lines %>% filter(name.fr == "La Seine")) %>%
#  st_transform(st_crs(iris.75))

# Export road and river layers to shapefile
st_write(roads.geom, dsn = "data/osmdata/road.shp")
st_write(river.geom, dsn = "lecture/data/osmdata/river.shp")
# bbox est utilisé pour centrer sur Paris
bb <- st_bbox(iris.75)
par(mar=c(0.2,0.2,1.4,0.2), bg="ivory")
plot(st_geometry(iris.75), col = "ivory", border="ivory3", 
     xlim = bb[c(1,3)], ylim =  bb[c(2,4)])
plot(st_geometry(roads.geom),col="#666666",lwd = 1.2,add=TRUE)
plot(st_geometry(river.geom),col="#87cdde",lwd = 3,add=TRUE)
plot(accidents.2019.paris %>% filter(grav==3) , pch = 20, col = "orange", add=TRUE, cex = 1)
## Warning in plot.sf(accidents.2019.paris %>% filter(grav == 3), pch = 20, :
## ignoring all but the first attribute
plot(accidents.2019.paris %>% filter(grav==2) , pch = 20, col = "red", add=TRUE, cex = 1)
## Warning in plot.sf(accidents.2019.paris %>% filter(grav == 2), pch = 20, :
## ignoring all but the first attribute

Géocodage

# avec banR
library(banR)
geo_banR = accidents.2019.paris %>% 
  filter(catv %in% c("VAE","EDP à moteur")) %>% slice(1:10) %>% 
  geocode_tbl(adresse = voie,code_insee = com) %>% select(latitude,longitude) %>%
  st_as_sf(coords = c("longitude", "latitude"),crs = 4326, agr = "constant") %>%
  st_transform(2154)
## Writing tempfile to.../tmp/RtmpyibIf9/file1fba837f8c41d.csv
## Warning: The `path` argument of `write_csv()` is deprecated as of readr 1.4.0.
## Please use the `file` argument instead.
## If file is larger than 8 MB, it must be splitted
## Size is : 666 bytes
## SuccessOKSuccess: (200) OK
## New names:
## * geometry -> geometry...10
## * geometry -> geometry...13
library(tidygeocoder)
## 
## Attachement du package : 'tidygeocoder'
## Les objets suivants sont masqués depuis 'package:banR':
## 
##     geocode, reverse_geocode
geo_tidygeocoder = accidents.2019.paris %>% 
  filter(catv %in% c("VAE","EDP à moteur")) %>% slice(1:10) %>%
  mutate(addr = paste(voie,", Paris, France")) %>%  
  geocode(addr,method="osm") %>% select(lat,long) %>%
  st_as_sf(coords = c("long", "lat"),crs = 4326, agr = "constant") %>%
  st_transform(2154)

st_distance(geo_banR,geo_tidygeocoder,by_element = TRUE)
## Units: [m]
##  [1] 950.52141 382.66079 122.00239 181.62935  47.91376 642.02833  41.90755
##  [8] 305.16959 664.76010 664.76010

Faire des cartes interactive avec R

De nombreuses solutions existent pour faire des cartes avec R :

  • ggplot2 est un package très utilisé pour faire tous types de graphiques, et a été adapté spécifiquement aux cartes (geom_sf).
  • Le package tmap contient des fonctionnalités avancées basées sur la logique de ggplot2
  • map_sf (anciennement cartography) s’appuie sur un langage dit “base R” et permet de faire des représentations cartographiques, basiques comme avancées.
  • mapview, leaflet et mapdeck permettent de faire des cartes interactives.

Par simplicité, nous nous concentrons ici sur ggplot2 pour la partie statique et mapview pour la partie interactive.

Petite introduction de sémiologie graphique

Light Semiology

Cartes interactives mapview

Les cartes interactives ne sont pas forcément très pertinentes pour représenter des informations géostatistiques.

En revanche, elles sont utiles pour explorer les bases de données. Voyons un exemple avec mapview concernant les accidents mortels à Paris en 2019.

#remotes::install_github("r-spatial/mapview")
library(mapview)
mapviewOptions(fgb=FALSE)
# construire une carte avec certaines options pour les cercles
# avec mapview la taille des cercles reste constante quel que soit le zoom. 
# grav = 2 : individus tués
mapview(accidents.2019.paris %>%
          filter(grav==2))

On customise un peu…

mapview(accidents.2019.paris %>%
          filter(grav==2),
        map.types = "Stamen.TonerLite", legend=FALSE,
        cex=5, col.regions="#217844", lwd=0, alpha=0.9)

On customise encore un peu plus…

mapview(accidents.2019.paris %>%
          filter(grav==2) %>%
          mutate(age=2019-an_nais),
        map.types = "Stamen.TonerLite", legend=FALSE,
        cex="age", zcol="sexe", lwd=0, alpha=0.9)

Cartes statiques ggplot2

ggplot2

Cartes avec ronds proportionnels

library(ggplot2)
ggplot() +
  geom_sf(data = iris.75,colour = "ivory3",fill = "ivory") +
  geom_sf(data = river.geom, colour = "azure",size=2) +
  geom_sf(data = roads.geom, colour = "#666666",size=0.5) +
  geom_sf(data = iris.75 %>%  st_centroid(),
          aes(size= nbacc), colour="#E84923CC", show.legend = 'point') +
  scale_size(name = "Nombre d'accidents",
             breaks = c(1,10,100,200),
             range = c(0,5)) +
  coord_sf(crs = 2154, datum = NA,
           xlim = st_bbox(iris.75)[c(1,3)],
           ylim = st_bbox(iris.75)[c(2,4)]) +
  theme_minimal() +
  theme(panel.background = element_rect(fill = "azure",color=NA)) +
  labs(title = "Nombre d'accidents de la route à Paris par iris",
       caption = "fichier BAAC 2019, ONISR\nantuki & comeetie, 2021",x="",y="")
## Warning in st_centroid.sf(.): st_centroid assumes attributes are constant over
## geometries of x

Cartes avec choroplethes

library(RColorBrewer)
bks <- c(0,round(quantile((100*iris.75$nbacc_blessgravtues/iris.75$nbacc)[iris.75$nbacc_blessgravtues!=0], na.rm=TRUE, probs=seq(0,1,0.2)),0))
pal <- c("#FFFFFF",brewer.pal(length(bks)-1,"Reds"))

ggplot() +
  geom_sf(data = iris.75,colour = "ivory3",fill = "ivory") +
  geom_sf(data = iris.75, aes(fill = 100*nbacc_blessgravtues/nbacc), colour = "grey80") +
  geom_sf(data = river.geom, colour = "azure",size=2) +
  geom_sf(data = roads.geom, colour = "#666666",size=0.5) +
   scale_fill_gradientn(name = "Part (En %)", 
                       values=bks/max(bks),
                       colours = pal) +
  coord_sf(crs = 2154, datum = NA,
           xlim = st_bbox(iris.75)[c(1,3)],
           ylim = st_bbox(iris.75)[c(2,4)]) +
  theme_minimal() +
  theme(panel.background = element_rect(fill = "azure",color=NA)) +
  labs(title = "Part des Accidents avec blessés graves et morts par iris à Paris",
       caption = "fichier BAAC 2019, ONISR\nantuki & comeetie, 2021", x="",y="")

library(RColorBrewer)

acc = iris.75 %>% 
  st_join(accidents.2019.paris) %>% 
  group_by(CODE_IRIS) %>% 
  summarize(nb_acc=n(),nb_vl=sum(if_else(catv=="VL seul",1,0)))


bks <- round(quantile(100*acc$nb_vl/acc$nb_acc, na.rm=TRUE, probs=seq(0,1,0.2)))
bks
##   0%  20%  40%  60%  80% 100% 
##    0   19   33   46   57  100
acc = acc %>% mutate(txaccvl = cut(100*nb_vl/nb_acc,bks)) 
pal <- c("#FFFFFF",brewer.pal(length(bks)-1,"Reds"))

ggplot() +
  geom_sf(data = iris.75,colour = "ivory3",fill = "ivory") +
  geom_sf(data = acc, aes(fill = txaccvl)) +
  geom_sf(data = river.geom, colour = "#87cdde",size=2) +
  geom_sf(data = roads.geom, colour = "#666666",size=0.5) +
  scale_fill_brewer(name = "Part (En %)",palette = "Reds", na.value = "grey80") +
  coord_sf(crs = 2154, datum = NA,
           xlim = st_bbox(iris.75)[c(1,3)],
           ylim = st_bbox(iris.75)[c(2,4)]) +
  theme_minimal() +
  theme(panel.background = element_rect(fill = "ivory",color=NA)) +
  labs(title = "Part des Accidentés velos par iris à Paris",
       caption = "fichier BAAC 2019, ONISR\nantuki & comeetie, 2021", x="",y="")

catv_ol = accidents.2019.paris %>% st_drop_geometry %>% count(catv) %>% arrange(n) %>% pull(catv)
gg = accidents.2019.paris %>% mutate(catv_o = factor(catv,levels=catv_ol)) %>% filter(catv_o %in% tail(catv_ol,10))
ggplot()+geom_bar(data = gg ,aes(x=catv,group=sexe,fill=sexe))

ggplot()+geom_bar(data = gg,aes(y=catv_o,group=sexe,fill=sexe))+theme_bw()+scale_fill_brewer("Sexe",palette="Set1")+labs(title="Nombre d'accidentés par type de véhicule et sexe", subtitle="à Paris en 2019, pour les hommes et les femmes ",caption = "fichier BAAC 2019, ONISR\nantuki & comeetie, 2021", x="",y="")

Crédits et reproductibilité

Présentation faite grâce au package rmdformats.

Elle s’inspire très fortement, ainsi que son tutoriel, d’une précédente formation donnée par les mêmes auteurs avec Timothée Giraud.

Partage de la configuration de R et des packages utilisés :

sessionInfo()
## R version 4.1.0 (2021-05-18)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
## 
## locale:
##  [1] LC_CTYPE=fr_FR.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=fr_FR.UTF-8        LC_COLLATE=fr_FR.UTF-8    
##  [5] LC_MONETARY=fr_FR.UTF-8    LC_MESSAGES=fr_FR.UTF-8   
##  [7] LC_PAPER=fr_FR.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=fr_FR.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] RColorBrewer_1.1-2 ggplot2_3.3.3      mapview_2.9.9      dplyr_1.0.6       
## [5] sf_0.9-8          
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.6              svglite_2.0.0           lattice_0.20-44        
##  [4] png_0.1-7               class_7.3-19            leaflet.providers_1.9.0
##  [7] assertthat_0.2.1        digest_0.6.27           utf8_1.2.1             
## [10] R6_2.5.0                leafpop_0.1.0           stats4_4.1.0           
## [13] evaluate_0.14           e1071_1.7-7             highr_0.9              
## [16] pillar_1.6.1            rlang_0.4.11            uuid_0.1-4             
## [19] raster_3.4-10           rmarkdown_2.8           labeling_0.4.2         
## [22] webshot_0.5.2           stringr_1.4.0           htmlwidgets_1.5.3      
## [25] munsell_0.5.0           proxy_0.4-25            compiler_4.1.0         
## [28] xfun_0.23               pkgconfig_2.0.3         systemfonts_1.0.2      
## [31] base64enc_0.1-3         htmltools_0.5.1.1       tidyselect_1.1.1       
## [34] tibble_3.1.2            bookdown_0.22           codetools_0.2-18       
## [37] fansi_0.5.0             crayon_1.4.1            withr_2.4.2            
## [40] grid_4.1.0              jsonlite_1.7.2          satellite_1.0.2        
## [43] gtable_0.3.0            lifecycle_1.0.0         DBI_1.1.1              
## [46] magrittr_2.0.1          units_0.7-1             scales_1.1.1           
## [49] rmdformats_1.0.2        KernSmooth_2.23-20      stringi_1.6.2          
## [52] farver_2.1.0            leaflet_2.0.4.1         sp_1.4-5               
## [55] ellipsis_0.3.2          brew_1.0-6              generics_0.1.0         
## [58] vctrs_0.3.8             tools_4.1.0             leafem_0.1.6           
## [61] glue_1.4.2              purrr_0.3.4             crosstalk_1.1.1        
## [64] yaml_2.2.1              colorspace_2.0-1        classInt_0.4-3         
## [67] knitr_1.33